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Abstract: We show that the solutions of the Bogomolny equations for the Abelian Higgs 
model on a two-dimensional torus, can be expanded in powers of a quantity e measuring 
the departure of the area from the critical area. This allows a precise determination of the 
shape of the solutions for all magnetic fluxes and arbitrary position of the Higgs field zeroes. 
The expansion is carried out to 51 orders for a couple of representative cases, including the 
unit flux case. We analyse the behaviour of the expansion in the limit of large areas, in 
which case the solutions approach those on the plane. Our results suggest convergence all 
the way up to infinite area. 
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1. Introduction 

Topological defects play an important role in Particle Physics, Cosmology, and many areas 
of Condensed Matter Physics, like superconductivity, superliquid helium, etc. Minimum 
energy(action) configurations carrying topological charges arise as solutions of non-linear 
partial differential equations. In some situations, these solutions are not explicitly known, 
and one needs to make use of approximate analytical or numerical methods to study their 
structure and properties. 

The Abelian Higgs model is one of the simplest models to study these ideas. It serves 
as a relativistic field theory extension of the Ginzburg-Landau description of a supercon- 
ductor, in which a scalar field represents the condensate of Cooper pairs. In addition to 
superconductivity, modifications of the Abelian Higgs model have found applications in 
other domains of Physics, such as Cosmology (see for example Ref. (1)). On a different 
level, the Abelian Higgs model acts as a simplified model in which to explore ideas and 
develop methods useful for non-abelian gauge theories. It illustrates phenomena such as 
spontaneously broken gauge symmetry, colour confinement, the dual Meissner effect, topo- 
logical charges, and others, all of which play a role in our present understanding of Particle 
Physics. 

Abrikosov (2) realised that superconductors contain string-like topological defects or 
configurations, which represent magnetic fiux tubes. A corresponding stable stationary 
cylindrically-symmetric configuration of minimum energy per unit length of the Abelian 
Higgs model was discovered by Nielsen and Olesen (3). Cutting a slice orthogonal to the 
direction of the string, we obtain a two-dimensional field configuration. Its stability arises 
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from the fact that it possesses non-trivial topological charge. In this case, the topological 
charge is just the magnetic flux going through the plane which, for finite energy (per 
unit length) configurations, is quantised in multiple units of lir. Despite its conceptual 
simplicity, there is no explicit analytical expression for the unit magnetic-flux minimum 
energy configuration (the Nielsen- Olesen vortex). 

The self-coupling constant of the scalar field A determines the ratio of photon to scalar 
field masses. There is a critical value of this coupling Ac (equal to one, in our units), 
separating the case of type I (A < Ac) and type II (A > Ac) superconductors (4; 5). 
Similarly to what happens with non-abelian gauge theories in 4 dimensions, Bogomolny (6) 
discovered that (for A > Ac) the 2-dimensional energy of any configuration is bounded by a 
given constant times the absolute value of the topological charge. At the critical coupling 
Ac, the bound is saturated by the solutions of a system of first order partial differential 
equations: the Bogomolny equations. 

Solutions of the Bogomolny equations not only exist but, for A = Ac, exhaust all 
solutions of the two-dimensional field equations (energy extrema) (7). Furthermore, the 
space of solutions defines a manifold of dimension 2q (8), where q is the number of flux 
quanta. This is again very similar to the situation for self-duality equations in non-abelian 
gauge theories in 4 dimensions (4D). Indeed, a connection between both topics is known (7). 
Despite these properties, there arc no explicit analytical expressions for these solutions. 
The unit flux cylindrically-symmetric solution was constructed as a power series in the 
radial coordinate by H. J. de Vega and F. A. Schaposnik (9). For higher fluxes, a solution 
is determined uniquely by the location of the q zeroes (counted with multiplicity) of the 
Higgs field. These points can be interpreted as the position in the plane of q unit- vortices. 
Since the energy does not depend on these positions (is given by the flux alone), one can 
think of these unit- vortices as non-interacting (This is not exactly true when one takes into 
account the kinetic terms in 4D). There are corresponding zero- modes associated with this 
degeneracy (10). This contrasts with the situation for non-critical coupling. Numerical 
simulations (11) and other methods tell us that when A > Ac, the interaction between 
vortices is repulsive, and there arc no static solutions of the equations of motion for flux 
greater that 27r in the plane. On the other hand, if A < Ac the interaction between vortices 
is attractive, and the solution of the equations of motion is a single vortex that carries all 
the flux (also called a "giant" vortex). 

Certain problems require some knowledge of the structure of the solutions, and, in 
the absence of analytical expressions for them, one must rely upon approximate methods. 
Frequently these methods are specific to cylindrically symmetric solutions, or make use of 
specific ansatze. This is the case of the numerical methods of Ref. (11), for example. Other 
approximations are based on expansions in powers of certain quantities. We already cited 
the study of de Vega and Schaposnik (9) for cylindrically-symmetric solutions. Similar ex- 
pansions have been used to study excitations of the vortices in the type II superconducting 
phase (12). 

In this paper we will present a new method to obtain the solutions of the Bogomolny 
equations on the 2-dimensional torus by means of an expansion in a parameter measuring 
the departure of the area from the critical area of 4g7r (measured in units of the square 
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correlation length). The method can be applied to obtain solutions of arbitrary flux and 
with Higgs zeroes at prescribed (but arbitrary) points. Truncating this expansion at a 
given order one obtains approximations to the solutions which tend to become worse with 
increasing area. However, as we will see, accessible truncations of the expansion give 
precise results even for large sizes, where the torus solutions approach those of the plane. 
Thus, although our primary goal is that of obtaining the solutions on the torus, it seems 
possible that problems pertaining to the plane can be addressed in this way too. Part of 
our motivation arises from the fact that a similar expansion was proposed to study self- 
dual Yang-Mills field configurations on the four dimensional torus (13). The present case 
is a simplified version of the 4D problem, and hence better suited for performing a more 
detailed convergence analysis. 

The outline of the paper is as follows. In the next section the method will be presented. 
In the following one, we will apply it to the two simplest cases: The unit-flux vortex case 
and a flux=2 case with no circular symmetry. In both cases the expansion is carried out up 
to 51 orders. This enables us to do an analysis of convergence of the series for various sizes, 
including the infinite area case. Finally, in the last section, we present the conclusions 
and indicate prospects for future research. In addition, we have included an appendix 
describing a variant method with a quantum mechanical flavour. Although, less efficient 
than the method explained in section 2, it has several advantages over it, one being that, 
at the present stage, seems better suited for generalisation to the non-abelian self-duality 
study. It also has served as a test of the actual computed coefficients of the expansion, 
since they have been obtained with both methods and two independent codes. 

2. Description of the method 

In suitable units (unit charge and photon mass), the Lagrangian density of the abelian 
Higgs model in four-dimensional(4Z)) Minkowski space-time, is given by: 



where (/> is a complex scalar field, and = d^i — is the covariant derivative with 
respect to the U{1) gauge potential. This system is known (2; 3; 8) to possess static, 
z-independent (vortex like) solutions of the classical equations of motion. These configu- 
rations are local minima of the 2D energy, whose density is 



where Latin indices run over two spatial coordinates: = 1, 2. 

We will focus in the case in which the coupling A takes the critical value Ac = 1. Then, 
the second order (2D) differential equations of motion reduce to a set of first order ones 
(Bogomolny equations) (6). For positive flux and in our units these equations take the form: 



(2.1) 



£ = ^^F,^F'i + \\D,<t>\^ + - If 



(2.2) 



(£>i +zL>2)<^ = 



(2.3) 
(2.4) 



B = -{l-\4>\'') 
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where B is the (z-component of the) magnetic field. 

Our goal is to study the solution to these Bogomolny equations for fields living in 
a 2D Torus (for an introduction to gauge fields on the torus see(14)). The appropriate 
mathematical description involves sections of a U(l) bundle. We will be working within a 
fixed trivialisation. The torus can be viewed as the quotient space of the plane M? modulo 
the lattice A generated by two linearly independent 2D vectors e*-^-*, e*-^-*. We assume that 
the torus is equipped with a flat Riemannian metric, which we will fix to be Euclidean. 
Within a specific trivialisation, the charged fields (sections) (j){x) = (l){xi,X2) are given by 
complex valued functions satisfying certain periodicity properties: 

(j){x + e'^^) = Qi{x)(l){x) (2.5) 

The U{1) fields r2j(x) are the transition functions, which have to satisfy the following 
consistency conditions 

ni{x + e(2)) Q2{x) = n2{x + e(^)) Qi{x) (2.6) 

The topological properties of the bundle, encoded in the transition functions, are associated 
with the first Chern class ci. Its corresponding integer Chern number is known, and will 
be shown later, to physically correspond to the magnetic flux going through the torus in 
units of 2tt. 

Without loss of generality we can choose the following specific form of the transition 
functions: 

Qiix) =exp{nTUj{e'^'\x)} (2.7) 

where to is an antisymmetric form. The consistency condition ( |2.6| ) forces 

uj{e^^\ e^'^^) = q to take integer values. This is precisely the first Chern number mentioned 

previously. 

Gauge fields are connections on this bundle. It is well known that the space of gauge 
fields is an affine space. The associated vector space is the space of 1-forms on the torus. 
Thus, we can decompose 

A = + a (2.8) 

where cj is a 1-form and A^^^ is a specific connection. For the latter, it is natural to select 
a gauge field having constant field strength / = F12 = —F21: 

Af\x) = -^F,,xj (2.9) 

Compatibility with the boundary conditions relates the antisymmetric matrix F with u) as 
follows 

F = 2TTUJ (2.10) 

which shows the relation between flux and Chern number. For the 1-form a we might use 
Hodge theorem and split it as a sum of of an exact, co-exact and a harmonic form: 

A = 7ruj{x,dx) +v + dg - 6h (2.11) 
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In components we might write 

Ai{x) = Af\x) + ^^g - eijdjh + (2.12) 

where h and g are real periodic functions and Vi are constants (a representative of the 
harmonic forms on the torus). Although locally these constants are pure gauges, globally 
they are not. Their value influences Polyakov loops winding around the torus, which are 
gauge invariant quantities. In this way it is easy to realise that the Vi can be considered 
elements of the dual torus M^/A*. In addition, one can fix h so that its integral over the 
torus vanishes (/ cFxh = 0). 

Now we can return to the Bogomolny equations and express them in terms of vi and 
the periodic function h. For simplicity we will restrict to orthogonal periods e^^^ = (/i,0), 
e(2) = (0, 12) (in appendix ^ we will study the general case). It is convenient to make use 
of complex coordinates: 

Xi + iX2 _ Xi-iX2 /o 10 \ 

z = ; z = (2.13a) 

d = d^-id2; d = di+id2 (2.13b) 

di = \{d + d); d2 = Ud-d) (2.13c) 
2 li 

The notation is chosen so that dz = dz = 1. 

Similarly the vector potential can be expressed as a complex function A = Ai — iA2 
(with A = Ai + iA2)- In this notation the Bogomolny equations take the form: 

(d-i'A)(^ = Q (2.14a) 

-'-[dA-dA) = \{l-\ct>\'') (2.14b) 

Our specific parametrisation of the vector potential (A*-*^-* = —if'z) becomes: 

A = dg-idh + v-ifz (2.15) 

where h and g are periodic functions, and v = vi — iv2 is a complex constant. For the 
Higgs field we will use the following parametrisation 

= Me~^+'3^ (2.16) 

where x is a normalised function (i.e. \x\'^d^x = lil2)-, satisfying the same boundary 
conditions as and N \s a. normalisation constant. With this terminology the Bogomolny 
equations become: 

{d + fz)x = wx (2.17a) 
ddh = \{l-2f-\M\^e~^''\x\^) (2.17b) 

These equations can be solved sequentially. The first equation allows one to determine x, 
satisfying the correct boundary and normalisation conditions. This can be done analytically 
as will be shown in the next subsection. Once this is solved, we can use the second equation 
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to solve for the periodic function h and the normahsation J\f. The equation is, however, 
a non-hnear partial differential equation and the analytic solution is not known. For a 
particular value of the area /1/2 = ^gvr, corresponding to / = ^, a solution is given by 
h = J\f = 0. Our strategy consists in using e = 1 — 2/ as a perturbative parameter. In 
what follows we will explain more precisely the procedure that we follow. 
Since h is periodic, we can use its Fourier decomposition: 

h=Y^ (2-18) 

nin2 

and expand the Fourier coefficients in a power series in e: 

00 

hn,n, = J2hi%,e' (2.19) 
fe=l 

Similarly, we can write 

00 

|AA|2 = ^Afce'= (2.20) 

k=l 

The only additional input that we will need are the Fourier coefficients of |xP- Although 
in some cases it is interesting and possible to deal with Fourier coefficients which are series 
expansions in powers e, we will restrict here to the case in which these coefficients are 
independent of e (but dependent on the aspect ratio r = The different situations 
will be clarified in the following subsection. Notice that our condition on h fixes /iqo = 0. 
Alternatively, we might have taken M = 1, and allow for non-zero hoo, but our choice is 
more natural within our expansion. 

To solve Eq. 2.17b| we equate the Fourier coefficients of both sides of the equation order 



by order in e. The left-hand side takes the form 

-il-enr,m,n2)Y.hi%^e' (2.21) 



where 



TTT 

^(T,ni,n2) = — 



„2 _L !!2 



(2.22) 



which vanishes for ni = 712 = 0. We remind the reader that r = ^2/^1 and q is the first 
Chern number. 

To treat the right-hand side we first expand it in powers of e. To order e the coefficient 
is given by (1 — Ai|xP)/2. The coefficient of order (for > 1) is given by 



A:=0 



N-k 



-2)* 



E 



/,(l)*^../,(iv-l)* 



Ixl 



(2.23) 
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Now, using the Fourier coefficients of |xp (see the fohowing subsection), we can obtain 
the Fourier coefficients of the expression above by applying a series of convolutions. We will 
label these Fourier modes with the symbol J-lvin2- These are functions of Aj(z = 1, . . . , N), 
and h^'''\i = 1, . . . , — 1). In order for the equation to have a solution it is necessary 
that J^iJ"^ = 0. This can be regarded as an equation for A^r. It allows to determine this 
coefficient uniquely in terms of Aj(z = 1, . . . , — 1), and /i*^*^(i = 1, . . . , — 1) 

Finally, the equation allows one to obtain the Fourier coefficients h^^^ to order as 
follows 



[ rii = 712 = 

This equation defines a recurrence which, starting from hn}n2 = 0, enables the determi- 
nation of the coefficients h^n}n2 uniquely. To order e, Tn^n2 are the Fourier coefficients of 
(1 — AilxP)- The vanishing of the ni = n2 = coefficient and the normalisation condition 
for X implies Ai=l. The coefficients hn}n2 ^"^^ then given by 



Computing the coefficients to higher orders demands performing convolutions, which 
involve infinite sums over several integers. We do not have closed analytical expressions 
for them. In practice, however, the Fourier coefficients decrease very fast with the order, 
so that a truncation of these sums to a finite subset allows, as we will show later, the 
numerical determination of the coefficients to machine precision up to high orders in the 
expansion. 

We emphasise that the previous procedure gives rise to a unique solution h. All the 
degrees of freedom associated to the space of solutions of the Bogomolny equations resides 
in the function x which will be treated in the following subsection. 



2.1 The function x 



In this section we will solve the first Bogomolny equation, Eq. 2.17a. First of all we will 
analyse the dependence on v = vi + iV2- It is easy to see that if Xv is the solution the 
equation for a value of v, then a solution for v' = v + 2ifa is given by 

X^, =ef^'''-''''^Xviz + a,z + a) (2.26) 

Since the prefactor is simply a phase, it is clear that solutions corresponding to different 
values of v represent solutions translated in space. Having this point in mind we can simply 
restrict from now on to = 0. 

To solve Eq. p. 17a we define rj = e^^'^^~^^x- Iii terms of this function the first Bogo- 



molny equation is equivalent to the condition of holomorphicity drj = 0. So rj = r][z) is 
an analytic function of the variable z. The boundary conditions for this function can be 
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obtained from the ones of (j) (Eq. ^.S] ), and are given by 



7] \ Z + l- 



1] [ Z + 



.k 



(2.27a) 



(2.27b) 



These are the typical conditions satisfied by theta functions (15). We will now analyse the 
space of solutions in various cases. 



For minimal flux q = 1 f = ^j^, it is easy to see that our function is given by 



2I2 
h 



(2.28) 



This, together with the normalisation condition, allows us to obtain the Fourier decompo- 
sition of x^- 



\X\ 



g ^ — ^ — — ^2mn-iXi/li ^2mn2X2/l2 |^_]^^'^l»^2 



(2.29) 



711,712 



where ^(t, ni,n2) is given by Eq. 2.22. From here it is trivial to obtain J-nin2 to initiate 
the iteration that determines h. 

For flux 27rq the space of solutions is multiple dimensional. There are several alternative 
ways to characterise individual solutions. One possibility is to fix the position of the zeroes 
of X) which coincide with those of the Higgs field (j). Then we can write 

f 2tt{z + Zc- UJl) ..l2\ a f 27r{z + Zc- UJ2) ..l2\ a f^T^iz + Zc - UJg) I2 
^ - [ k ' t J [ h K) [ h K 

(2.30) 

where Wj are complex constants, and z^ = \ih + ih)- This is a holomorphic function, and 
satisfies the correct boundary conditions provided 



^ = qzc 



(2.31) 



The function -q defined in this way has q zeros located at the points Wj. Eq. 2.31 then 
specifies that the centre of mass of all the zeroes is located at the centre of the torus 
(y, y). To shift the centre of mass one must choose f 7^ 0. 

An alternative description of the space of solutions follows naturally from the quantum 
mechanical formulation of appendix A basis of the space of holomorphic functions 
satisfying the same boundary conditions as 'q{z) is given by: 



r]s{z) = -d 



s/q 




2TTqz I2 



(s = 0,...,g-l) 



(2.32) 



where the symbols 



denote Theta functions with rational characteristics: 



{n+af ^2i{n+a)(z+b) 



(2.33) 
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The functions ris{z) have q zeros located at 



{k = 0,...,q-l) 



(2.34) 



An arbitrary solution of the problem is given by a suitably normalised linear combination 
of r]s{z): ijo^^s CsTlsiz). 

Now we can relate the two descriptions by relating the coefficients Cg to the position 
of the zeroes Wj. We can construct a linear combination that possesses q — 1 zeros at the 
(distinct) complex points LOi, i = 1, ... ,q — 1 hy setting 



ri{z) oc 



r]o{z) r]i{z) 



riq-i{z) 



(2.35) 



General theorems about elliptic functions tell us that the remaining zero of this r] function 
is located at a point ujq which enforces the centre-of-mass condition. 
The coefficients 

Cs = £sii...ig-i'nii{^l) ■ ■ ■ Vig-ii^q-l) (2.36) 

give rise to coordinates in the manifold of solutions of the Bogomolny equations. Actually, 
they define coordinates in the submanifold defined by the centre of mass condition for the 



zeroes Eq. 2.31. Using properties of elliptic functions and the result of Taubes (8) one 
can easily see that this submanifold is diffeomorphic to CP''"^ and Cs are homogeneous 
coordinates. One has to include v to give coordinates over the whole manifold of solutions. 
Since the latter varies over a torus, we recover the fiber bundle structure described in 
Refs. (16; 17; 18). 



The formula relating the coordinates to the zeroes Eq. 2.36 fails if two or more zeroes 



coincide. To write the correct formula one has to substitute some rows of Eq. 2.35 with 



the values of the derivatives (up to order given by the degeneracy) of the functions rjs . 

Finally, it is easy to obtain the Fourier coefficients of |xP, using the corresponding 
ones for theta functions with rational characteristics. The result is 



1 



\X\ 



2 e 1 



2m — ^ 

, e 1 



o ■ ^1 



ly^ 2mn2jr 
1 e '2 



(2.37) 



where Cs is given above. 

In all our previous construction the position of the zeroes of the Higgs field have been 
fixed relative to the Torus lengths /i, /2- One could, in principle, insist in fixing the position 
of the zeroes in units of the inverse photon mass. This amounts to taking the positions uJi 
as functions of e. We might still use the previous formula but now one has 



\X\ 



X e 



(2.38) 



The iterative solution of Eq. 2.17b given in the previous subsection must then be modified 
by combining the powers of e from |xp into the expansion. 
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3. Representative examples 

In this section we will show the effectiveness of the method presented in the previous 
section, by applying it to two explicit examples. The first one is the standard q = 1 case, 
which should converge towards the usual Nielsen-Olesen vortex with cylindrical symmetry 
when the Torus is big enough. An extensive analysis of convergence of the series for 
different volumes will be presented. The second example is a g = 2 case which has no 
remnant continuous spatial symmetry. This is important since many alternative methods 
are specific to cylindrical symmetry. 



We have applied the machinery explained in the previous section to the unit flux q = 1 
and unit aspect ratio t = I2/I1 = 1 case. The solution is unique up to translations. We 

(k) 

have obtained the coefficients of the expansion /inin2 up to order k = 51. Convolutions 
were performed by truncating the sums to Fourier modes in the range \ni\ < nmax- This 
was a implemented using a FORTRAN 90 program running in a standard PC. Results up to 

A; = 51 required computation times around 100 hours. 

We have, first of all, analysed the effect of truncation in the number nmax of Fourier 
modes used in convolutions on the value of the coefficients. We explored the cases nmax = 

10, 14 and 20. We estimate the relative difference (difference over sum) in the coefficient 

(k) 

hnin2 obtained from two different truncations 10-14 or 14-20. The difference is attributed 
to an error in the value for the smaller nmax- The differences increase with n, and with 
k. However, even for the highest order (k = 51) the relative difference between the results 
with nniax = 14 and 20 is of the order of machine precision 10^^^'"^'^ for all modes having 
Ui < 7. For higher modes it is to be expected that the error is mostly due to an error in 
the values for nmax = 14, so the modes for nmax = 20 remain probably within machine 
precision for higher n^. To analyse this we used the comparison 10-14 to try to understand 
the dependence of errors with n^, k and nmax- It is to be expected that the error is 
proportional to (some power) of the size of the neglected terms /limaxnmax- This fits nicely 
with the observed dependence of the relative difference 



Our interpretation of the source of the errors suggest that the coefficients multiplying | 
are proportional to nmax- This allows us to scale the results to nmax = 14- Indeed, our 
data on the difference 14-20 is consistent with this interpretation, although there are few 
values of k and n^ to allow an analysis by itself. On the basis of these results, we conclude 
that, even at A; = 51, the coefficients obtained for nmax = 20 are correct up to machine 
(double) precision for n^ < 14 — 15. 

Once the coefficients arc obtained, one can reconstruct the Fourier modes of h, the 
magnetic field B and the modulus of the Higgs field |0| for arbitrary torus areas A = I1I2 = 
Y^- Obviously, the precision of the truncated series depends on the value of e, decreasing as 
e increases up to the infinite area case of e = 1. Even for relatively large areas the shape of 



S.l q = l 




(3.1) 
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the reconstructed function looks quahtatively quite good. See for example Fig. 1,2, where 
we display the magnetic field for e = 0.9, obtained from the Fourier modes (|nj| < 20) 
computed with 51 orders in the expansion. 




10 20 30 40 50 



Figure 1: Shape of the magnetic field B{x) Figure 2: Contour plot of B{x). 

obtained for e = 0.9 with 51 orders in the 

expansion. 



To go beyond the qualitative level and estimate the accuracy of the truncated series, 
we use the degree of satisfaction of the Bogomolny equation as a measure of the error. 
Thus, we compute the magnetic field B{x) for various values of e and compare it with the 
right-hand side of Eq. ^(1 — |</)(x)p). The latter is computed using the truncated 

expansion of h and the parametrisation Eq. 2.16. More precisely we computed the L2 and 
Loo norm of the difference: 



L2{N,t) 




B(x) 



1 



1/2 



(3.2a) 



max 

X 



B{x)--il-\<p{x)\'] 



(3.2b) 



where is the maximal order in the expansion. We have analysed the and e dependence 
of both quantities. Results are qualitatively the same for both, so we will choose L2 to 
display. First, we will comment on the maximum precision, attained for n = 51. For e < 0.6 
the L2 norm is compatible with zero within machine precision (order 10"^^'^^''). Beyond 
this value L2 becomes sizable and increases, reaching 10""^ at e = 0.95. For comparison we 
point out that with = 1 the value (at e = 0.95) is O(10~^). Convergence is therefore 
slow in this case, but notice that the linear size of the box is 15 times the Debye screening 
length or 4.5 times the square root of the critical area. 

We performed a more systematic study by analysing the N dependence for fixed value 
of e, in the range 0.1 — 0.95. In this range the results are unaffected by the truncation 
in the number of Fourier modes rijnax- In all cases, we found that, beyond the first few 
orders, the dependence of L2{N,e) with N oscillates around an exponential fall-off. As an 
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Figure 3: L2{N, e) versus (in log^o scale) Figure 4: a(e) versus e compared to our best 
for e = 0.9. fit Eq. 3.5. 



example, we show in Fig. 3 the e = 0.9 case . Therefore, we fitted log^g -^2 data to the 
following linear function: 

logioi^2(iV,e) = a(e)iV + 6(e) (3.3) 

The parameter a(e) is determined with errors reflecting the statistical and systematic un- 
certainties (range of fitting for example). Its value determines how the approximation 
improves when increasing the order A'^ in the expansion. Its negative value is an indication 
that the expansion is indeed convergent. Obviously as e increases so does a(e). This is 
displayed in Fig. 4. For a convergent series and small e one expects 

logio L2{N, e)^{N + 1) logio(e) + logio (civ) (3-4) 

and hence, a{e) = log^g(e) + constant + 0(e). This is indeed the behaviour shown by the 
Fig. 3. Fitting o(e) to a linear function of log^o e gives: 

a(e) = 1.022(5) logio e - 0.0106(6) (3.5) 

Errors reflect the quality of the fit. Remarkably nothing seems to be happening at e = 1, 
where the area diverges. Data cannot be taken directly at e = 1 because they are severely 
affected by the truncation in the number of Fourier modes, but the behaviour up to e = 
0.95 shows no sign of a change of pattern and extrapolates to a(l) < 1. Similar smooth 
behaviour is shown by 6(e). So wc take our results as an indication that the series actually 
converges all the way up to e = 1. 

Within our spirit of identifying the L2 (or Loo) norm of the equation with the error 
on the Higgs and magnetic field, we can use our data from a more practical viewpoint as 
an estimate of the number of terms required in the expansion to attain an a priori decided 
precision. An approximate formula can be derived from our data. If one is willing to 
compute the magnetic field with an error of 10"^ then the number of terms required in the 
expansion is given by: 

_ 1.83e + p-3.25 

0.01 - 1.02 logio e ' 

Although the formula gives a finite number even for e = 1, we stress once more that in 
practise at that very large volumes, truncation in the number of Fourier of terms would 
make the expansion increasingly computationally costly. 
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Now we will explore the implications of our expansion for large volumes. Our main 
assumption is that the solutions on the torus do converge to those on the plane. The 
convergence is expected to be fairly fast. A torus configuration is equivalent to a periodic 
array of vortices on the plane. However, vortices are exponentially localised objects so that 
if the period is large compared to the typical size of a vortex, the effect of the replicas is 
presumably very small. Now the convergence of the solution implies the following behaviour 
of the Fourier modes: 

Y{P) ^ 1^ = lime(l,ni(p),n2(p))/i„,(,),„,(,) (-1)"^^)+-^^ (3.7) 

where B{p) is the Fourier transform of the magnetic field on the plane. The limit is taken 
at fixed p given by: 

p'^lpf = e(l-e) (3.8) 
where ^ is given by Eq. |2]2|. This means that as e tends to 1, the integers ni{p) have to 
grow. If instead, we take the limit e ^ 1 keeping rij fixed, the values should converge to 
y(0) = 0.5 irrespective of rij. In our case (g = 1), computing the value at e = 1 using our 
51 orders and rii = l,n2 = we get 1^(0) = 0.499947172199. Worse results follow for higher 
modes (0.499572 for m = n2 = 1, 0.465354 for m = 2, n2 = 1, 0.387814 for m = rig = 2, 
etc). The numerical agreement provides an additional hint that the expansion converges 
up to e = 1. It also indicates a poorer convergence for larger rii (see later). 

For p non-zero, Eq. 3^ and the expectation of fast convergence, suggests that tuning e 
and TLi in such a way that p is fixed we should obtain similar values. Only p, the modulus of 
p, matters due to the cylindrical symmetry of the q = 1 solution on the plane. This is also 
satisfied by our expansion to a fairly high precision. For example, Y{p) can be computed 
for p'^ = 7r/20 using e = 0.975 and ni = n2 = 1, or e = 0.95 and ni = 1, n2 = 0. From our 
expansion we get 0.39967 and 0.39977 respectively. This number is presumably very close 
(< 1%) to Y{p) on the plane. Similarly for / = vr/10 we get 0.326083 and 0.326075 from 
the same two modes and e = 0.95,0.9 respectively. For p^ = tt/5 we get 0.22627 (e = 0.8), 
0.22676 (e = 0.9) and 0.22608 (e = 0.95). In this way we can use our expansion to compute 
the Fourier transform of the magnetic field for a vortex on the plane with a precision of a 
few percent. 

Now we will try to extract the consequences of the convergence of the expansion for 
the Fourier modes to a universal function of p as e ^ 1. For very large rii (large .^) we can 
compute Y{p) by taking e = 1 — p'^/C- Thus, 

TV 

Y(p) ^ _^(_l)-l+«2/,^^^^(,) ^ _^^(_l)ni+n.;,(fc^g-P^| (3.9) 



A;=0 



This suggests that for ^ — > cxd 



^SU-^^^(^/o(-ir+"^ (3.10) 

where the function ip{y) satisfies 

/"OO 

Yip)= / dy^iy)e-P"y (3.11) 
Jo 
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Fourier transforming we get: 



dy 

y 



ip{y)e 



-\x\y{4y) 



(3.12) 



One can test these considerations by computing approximants to (f{y) by using Eq. 3.10 

(k) 

for finite n^. In Fig. 5 we show the shape obtained from the coefficients /inin2- For any y 
we plot only those values of rii such that ^ > y/25. The smoothness and small dispersion 
of values agrees with our expectations. We also investigated the way in which the limit 
is approached for large ^. For example in Fig. 6 we plot — ^^(— l)"^'''"'^/ii'^n2 for different 
values of and a fixed value of y = k/^ = l/vr. The solid line is the result of a fit to 
a function of the form a + |. Similar behaviour obtains for other y values. This analysis 
could be used to obtain a more precise estimation of the value ip{y). For the time being 
we simply used the non-extrapolated shape shown in Fig. 5 and analysed the behaviour 
for large and small values of the argument y. For small y, the function is well described 
by exp{a' — with a' and h' very close to 1. For large y the behaviour is also very well 
described by an exponential exp{— a — hy}. A fit in the range y G [2,6] gives a = 0.2973 
and h = 0.9443. Assuming our formula Eq. |3.12 , we can, by saddle point methods, relate 
the large |x| behaviour of B{x) to these parameters. Indeed, h is predicted to be 1. The 
parameter a is given by — log(Zi/2) where Zi was obtained numerically by de Vega and 
Schaposnik (Zi = 1.7079)(9), and recently Ref. (19) predicted its value to be log(2)/4. 
These values of Zi/2 differ by 10% from e^". This is a quite satisfactory agreement for 
the non-extrapolated curve obtained from the coefficients of our expansion. 



0.5 



1.5 2 2.5 3 3.5 



Figure 5: ip{y) vs. y computed using 
Eq. 3.10 



Figure 6: -C^(-l)"i+"^/il';U for different 
values of ^ and y — k/^ ^ I/tt fixed. 



From formula 3.12 one can deduce a connection between our expansion and that of 
Ref. (9) in powers of The expression becomes 



S! 



oo 1 



(3.13) 



Numerically integrating the data we get Di = -0.999976, D2 = 0.747034, = -0.523573, 
D4 = 0.36505, in good agreement with Ref.(9) (-1, 0.72791, -0.48527, 0.31444 respec- 
tively) . 

From all the discussion above we see that our expansion, though originating from a 
small volume expansion on the torus, matches nicely with results known for the single 
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vortex case at infinite area. Even though it might not give the same level of precision 
as other methods in that regime, it has the advantage of being readily generalisable to 
arbitrary fluxes and positions of the Higgs field zeroes. Furthermore, the same coefficients 
provide solutions in arbitrary torus sizes. 

We can actually employ the previous formulas to estimate the error committed in 
hn^n^{e), the Fourier coefficients of the function h, as a result of the truncation of the 
series. The contribution of terms higher than in the expansion, Ajv/i„^„2 can be 



estimated in terms of the function ^p{k) and Eq. 3.10. The appropriate formula is 

/51/5 



(-1)-" A„.„.„(,) = j iy . ^^^-^ (3.14) 



The last equality is obtained from the asymptotic behaviour of ^{y) and, thus, is only valid 
for (n^+rig) < 10. Applying this formula we get results which match with the discrepancies 
observed in some cases. For example, as we said before, 1^(0) (e = 1) should be equal to 0.5 



irrespective of which mode is used to compute it. However, our formula Eq. 3.14 predicts 
that the truncated evaluation up to 51 orders and ni = 77-2 = 2 should fall short by 0.1105. 
The actual discrepancy found previously is 0.1122. Everything fits nicely. Our formulas 
can also be used to estimate the number of terms in the expansion required to obtain a 
given Fourier coefficient on the torus with a certain precision. 

3.2 q = 2 

Here we will apply our method to a multivortex situation. We take unit aspect ratio 
(r = 1) and two units of flux. We can also use the procedure explained previously to fix 
the position of the zeroes of the Higgs field. We took the following points: 

0.35/1, |V f 0.65/1, |V (3.15) 



2y V 2, 

which are separated along the x direction a distance 0.3/i. 

In Figs. 7,8 we display the shape of the magnetic field obtained for e = 0.9 and 51 
orders in the expansion. There is no particular difference in computational cost between 
this case and the unit flux one. The effects of truncation are similar to those obtained in 
the previous section: modes up to max(ni, 712) < 15 are calculated up to machine precision. 
Note, however, that the position of the zeroes introduces a new scale in the problem, which 
translates into a typical scale for the modes. This might cause trouble if the zeroes are 
very close together. 

We have repeated all our previous analysis of convergence with qualitatively identical 
results. For example, the L2 norm of the function {B{x) - ^(1 - |(/)|2)), noted L2{N,e), 
seems to fall off exponentially fast with A^, as in the q = 1 case. With similar definitions 
and methods to the ones used for q = 1 we got a(e) and 6(e). Our best fit to the former 
quantity now gives: 

a(e) = 1.032(6) logio e - 0.020(1) (3.16) 
Our previous conclusions about convergence extend to this case as well. 
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10 20 30 40 50 

Figure 7: B{x) field for e — 0.9, and 51 Figure 8: Contour of B{x) for e — 0.9, and 
orders. 51 orders. 

Focus It is interesting to compare the precision of our results with those employing 
other methods. In particular, the shape of the solution was computed 

Focusing on the results for large areas, we emphasise that the main effect of a change 
in e will the be to vary the separation between the vortices. Thus, with the coefficients 
obtained from our analysis we can actually explore the variations in shape for nearby 
vortices as a function of separation, a study which can have some interest (see for example 
Ref.(20)). In this (e ~ 1) case, we can also compare with other alternative methods. In 
particular, in Ref. (21) the two vortex solution is computed by numerical methods. The 
finite square size used corresponds to e = 0.91 and the precision attained 10~^. Our 
formulas give a precision in the 10"^ - 10"^ for this size, which is, at least, as good. 

4. Conclusions 

Let us summarise our results. We have shown how one can expand the solutions of the 
Bogomolny equations on a two dimensional torus in powers of e = 1 — 2/, where / is the 
average magnetic field (flux over area). The coefficients of the Fourier modes can be con- 
structed using an iterative procedure involving convolutions. Although, no close analytical 
expression for the coefficients exist beyond the first non-trivial order, these coefficients can 
be determined up to double precision machine accuracy (15-17 significant digits), by trun- 
cating to a finite number of modes. This method can be applied to construct solutions with 
arbitrary flux and location of the Higgs field zeroes. We have computed the coefficients for 
a couple of cases {q = 1 and q = 2) and the results are very encouraging. The 51 order trun- 
cated expansion is estimated to describe the shape of the function within machine precision 
up to areas which are two and a half times the critical area. But meaningful values can be 
extracted also for large sizes, where due to the exponential localisation of the solution, the 
configurations are close to those of infinite area. In particular, these results match nicely 
with what is known about the unit-vortex on the plane. Turning the information around, 
this allows to obtain precise expectations about the behaviour of the coefficients for large 
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order, which are satisfied by the data. No significant difference in performance is observed 
when studying multivortex solutions. 

The solutions on the torus are relevant to depict the behaviour of the system in a 
situation of high vortex density. Their thermodynamics was analysed in Ref. (16). The 
description in terms of Fourier modes has been used previously in other contexts, like in 
the study of skyrmion crystals Ref. (22; 23). It seems, however, on the basis of our results, 
that our method can be used successfully to study the infinite area case as well. 

There are a number of possible applications and generalisations of the method to other 
problems or situations, some of which are currently under study. Special mention deserves 
the application to self-dual configurations on the four dimensional torus. In this case, there 
are no explicit analytic formulas for the solutions and the present method might provide 
good results, as an alternative to numerical methods (24). Actually, the first term in the 
expansion was obtained previously (13) and served as initial motivation of this work. 

There are other interesting problems for which the present method can be used. For 
example, in the study of the dynamics of vortices, specially in the low energy limit in which 
the geodesic approximation is valid (25; 26). The main issue here is the determination of the 
metric within the manifold of solutions. This metric can be extracted from the behaviour of 
the solution itself in the vicinity of the Higgs zeroes (21), which suggests that our method 
can be successfully used(27). This study offers the opportunity to express and analyse the 
conserved quantities studied in Ref. (28) within our formalism. 

Finally, it is interesting to notice that the critical has a generalisation to 

Higgs-gauge systems in any Kahler manifold, in what is called the Bradlow limit(29). From 
this point of view, our expansion can be viewed as a particular case of a much more general 
concept: an expansion in the Bradlow parameter. Although part of our technology is 
specific to the torus, we think that there are appropriate generalisations to other manifolds 
and Riemannian metrics by using a different set of basis functions. Indeed, the lowest term 
has already studied for the case of the two-sphere(30). These generalisations are currently 
under study by the present authors. 



A. Quantum mechanical formalism 

In this section we will describe the formalism in quantum mechanical terms. Here we 
follow the spirit, notation and formulas of Ref. (31). Our goal is to characterise the space 
of sections of a U(l) vector bundle on the two-dimensional torus within a fixed trivialisation. 

Fields satisfying the boundary conditions (^^2.7) make up a pre-Hilbert space Hq 
with scalar product 

($1$) = !/' dx <^*{x)'^{x) (A.l) 
A Jt^ 

where A is the area of the torus. The integration is over the torus and because of the 
integrand periodicity can be performed over any fundamental cell. 

Following the standard quantum mechanical procedure we will now look for a complete 
set of commuting operators which can serve to find a basis of the Hilbert space. One family 
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of operators Ua is labelled by elements a of the dual lattice A*: 

Ua*(x) = exp{27rza(x)}*(x) (A.2) 

where o is a linear function of x satisfying a(e(*)) G Z. Notice that we can introduce two 
special elements Ui of A* associated with the form lo: 

Uiix) = ^ ^' ' (A.3) 

The set of these two elements is the dual basis to {e^^^}. All the operators are mutually 
commuting. They arc just gauge transformations of a special kind. 

In addition we will look at operators implementing translations. However, ordinary 
translations map out of Hq because the transition functions depend on x. It is clear on the 
basis of the nature of our fields that we should replace translations by appropriate parallel 
transporters. For that purpose we will make use of the privileged connection A^^^ on the 
torus, having constant (or uniform) field strength: 

= 7ra;(x, dx) (-A-.4) 

Finite translation operators 7^ can be defined as parallel transporters along straight 
lines: 

(Ta*) = e"'^7 + a) = e-*^'"(^'«) ^{x + a) (A.5) 

Obviously these operators do not commute: their commutator is determined by the flux of 
the gauge field through the corresponding parallelogram. 

Our boundary conditions can be then reformulated by saying that elements of Hq, are 
those fields left invariant by the operators: 

'?'e(0U(_5^^) (A.6) 

We can then introduce the operators 

ir«=r,w/,U^, (A.7) 

and re-express the boundary conditions by saying that (K^^^^ = I. Furthermore the 
operators obey: 

if(i)i;f(2) = exp{27rz/g} K^'^^K^^^ (A.8) 

The operator K^^'> generates a Zq group and its eigenvalues are given by exp{— 27rz|}, where 
s is an integer modulo q. Thus, one can decompose Hq into the q orthogonal subspaces of 
eigenvectors: 

Hq = ®UHq,s (A.9) 

K(2) 

maps Hq^s into Hq^s-i- The translation operators 7^ commute with K^*) and hence 
leave these subspaces invariant. 

Our task of finding a complete set of operators is achieved by K^^^ and an operator 
involving translations. As we will see in the study of the Bogomolny equation it is natural 
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to select this operator to be the covariant Laplacian constructed with the constant field 
strength gauge potential. Up to now the choice of metric in space has played no basic 
role. Here, however, this operator depends on the metric. We will take the metric to be 
Euclidean. Within this metric the lattice vectors e^*-* have a well defined length and scalar 
product. It is always possible to make a coordinate transformation to bring e^^^ to the 
form (/i,0). The other vector e*-^^ is then given by l2{cosip,smip). This reduces for <p = 
to the special case considered in section ^. 

The generators of translations along each axis are precisely the components of the 
covariant derivative corresponding to the A^^^ field. We obtain 

=di+ TTiqeijXj/A (A. 10) 

where €ij is the antisymmetric tensor with two indices(ei2 = 1), and A is the area of the 
fundamental cell. These operators are anti-hermitian and satisfy the following commutation 
relations: 

|i,<»),o<»)] = -./.zHp (A.n) 

After an appropriate rescaling this is just the Heisenberg algebra satisfied by momentum 
and position operators in one dimensional Quantum Mechanics. Using standard formu- 
las one can construct operators a, a''^ satisfying the commutation relations of creation- 
annihilation operators, and express the covariant derivatives in terms of them: 

i?f) = iy|(a + at) (A.12) 

Z)f = y|(a-at) (A.13) 

The covariant Laplacian associated to the A^^^ field is proportional to the Hamiltonian of 
a harmonic oscillator: 

i?fi?f = -/(ata+i) (A.14) 

Thus the space of classical sections of a U(l) bundle, has identical structure as the Landau 
levels of the quantum system. 

Finally, a basis of our space of sections is provided by the states |n, s) which are 
simultaneous eigenstates of the number operator and K^^\ where n is an arbitrary non- 
negative integer and s an integer modulo q. We have the following relations 

K(i)|n,s) = e~'^|n,s) (A.15) 
K(2)|n,s) = |n,s - 1) (A.16) 

Df^\n,s) = i^{V^^^\n + l,s) + V^\n-l,s)) (A.17) 
|n, s) = {-V^^^^\n + l,s) + y/n\n - 1, s)) (A.18) 
We consider the \n, s) states to be orthonormal within the scalar product Eq. A.l 
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We will now give the explicit form of the functions ^n.sixi, X2) corresponding to these 
states (|n,s)). For that purpose notice that any function ^{xi,X2) G Tig can be expressed 
as: 

pGZ 

which follows from analysing the periodicity under e'-^^ Now imposing that ^ belongs to 
Tlq^s, one concludes that, in the sum appearing in Eq. A.19 , p is restricted to p = s mod q. 
Imposing now periodicity under e*-^-* we arrive at: 

^{xi,X2) = e-*^'?'^i(^)"2(^) e-2^*P'^2{^) J(quj^(x) + p) (A.20) 

p&s+qZ 

Now we look at eigenstates of the number operator. For that purpose it is interesting to 
look at the way in which the creation and annihilation operators act on the function J 
appearing in the previous formula. We have: 

(A.21) 



V2 ay' smy? 



{^\J){y) = -^{±-^^y')J{y) 
V2 ay' smv? 



(A.22) 



where y' 



2it12 sin y 

After a standard quantum mechanical calculation we arrive at the expression of the 
function ^ns{xi-,X2) corresponding to the state \n,s) 

pes+qZ 

where Hn is a Hermite polynomial and y' = \/~l(x2 + ) ■ To reduce to the case of 

orthogonal e*^*^ one can fix = ^, a;i(x) = X2//2 and ^2(2;) = —x\ll\. 
For the special case of the vacuum state (n = 0), we get 



h 



s/q 




iz,r) (A.24) 



where 'd 



s/q 




(z, r) is a theta function with rational characteristics with arguments: 



z = tt{tuji - quj2) 
r = e^-qf 



(A.25) 
(A.26) 



Now we need to deduce the action of the translation operator T(a) on these states. The 
translation operator is defined as T{a) = ex.p{aiDf'^ + 02.02''^} = exp{— (z(a)a — z*aL')}. 
Now we can compute the matrix elements 

{m,s'\T{a)\n,s) = J,,, e^^("-'^)(-l)(^'^+") x (A.27) 



e 2 \z\ 



\J n\m\ IzP-^ 



i!(M- j)!(j + |m-n|)! 
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where M = min (m, n) and 



z(a) = ^{ai - ia2) = Izle'f^ (A.28) 



We can also compute the matrix elements of the operators Ua- Expanding a in the dual 
basis of A* we can write d{x) = —kiui2{x) + k2U}i{x). We will then denote Uq = U{ki, ^2)- 
This operator can be expressed in terms of K^^") and translations as follows: 

C/(A;i,A;2) = e-'"^ (i^(i))fc2T(a(A:i, ^2)) (A.29) 

From here we can compute the matrix elements: 

^ j-,(M-j)!(j + |m-n|)! 
where a and ^ can be obtained from the complex number 

K = -i^{lik2 - l2kie'^)/q = \[i e"^ (A.31) 
and ^ is the generalisation of ^ defined in Eq. 2.22 to the case of non-orthogonal e^*^: 

i = \K\'^ = -^^{kj^-^ + '^hk2 cos(^) (A.32) 

q sin ip li I2 

A.l The Bogomolny equations 

We can now reformulate the Bogomolny equations in this formalism. We can write the 
Higgs field as = eip, where ■0 is a normalised element of 7iq. It can be decomposed in 
our basis 

00 

if) = '^^^c{n, s)\n, s) (A. 33) 

n=0 

where the sum of the modulus square of the coefficients c(n, s) equals unity. Similarly the 
function h appearing in Section ^ is expressed in terms of a generalised Fourier decompo- 
sition 

h{xi,X2) = /lfc^fc^e2"^(-'=l'^2(x)+tea.i(x)) (^_34) 

fcifc2 

The first Bogomolny equation can then be expressed as follows: 

aV' = i-^ + id2)h) ^ (A.35) 
V ^/ 



Hence, using the decompositions of h and ip and the matrix elements deduced in this 
Section, we can re-express the Bogomolny equations as follows: 

(1 - e)|K|^/ifeifc2 = I X] c*(m, s')c(n, s) (m, s'|C/(-/ci,-/c2)|n,s) (A.36) 

m,s' ,n,s 

y/n + lc(n + = - ^ K*hk^k2 c{p,s) {n, s'\U{ki, k2)\p, s) (A.37) 

kik2 P,s 
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where k = ^{kil2e~^'^ — k2h)/q. The volume dependence of these equations is exphcit 
(contained in the dependence on e). Our method consists then in expanding the unknown 
coefficients in powers of e: 

oo 

c{p,s)=Y,c^'\p,s)e' (A.38) 

fc=0 

oo 

hk.k, = J2€i'' (A.39) 

k=l 



and solving the equations order by order. 

Let us compute the leading terms for the q = 1 case. To lowest non-trivial order we 
have 

c(°^H = 5„o (A.40) 
and h^^^i^^ = 0. Then, plugging this value into the first Bogomolny equation we get: 

= -~m{-h - = e-«> (A.41) 



The following order obtains from the other equation 

1 1 ( K*^'"' 

cW(n) = --7^ E '^*h('Hh,k2) (n - l\U{h,k2m = --^ ^ -J- ^"^ (^-^2) 

kik2 k\k2 



This is real and vanishes for odd n. 



One can then continue to iterate Eqs. A.3(:- A.37 to higher orders in e. To evaluate 



the expressions one has to restrict to a finite number of elements on the basis. Now, in 
addition to the cut in Fourier modes n^ax) one has to cut in the spectrum of the number 
operator n < np^rt- As the order k grows both n^^^ and np^^t should grow to keep the 
numbers within machine precision. We have developed an independent FortranQO program 
to evaluate the c^^\n), h^^\ni,n2) using this procedure. With this program we have 
reproduced, within machine precision, the coefficients in the expansion of the low lying 
modes rii up to 30th order in the expansion for the q = 1 case. This serves as a check that 
there are no unexpected bugs in the determination of the coefficients. Increasing the order 
one starts noticing sizable errors associated to the cut in n^^^^. Unfortunately, increasing 
this number is not only limited by computer resources, but also by numerical instabilities. 
For example the computation of the matrix elements (n, s'\U{ki, k2)\p^ s) becomes unstable 
for large fej, n and p. This is due to wild cancellations of large numbers appearing in their 
definition. Actually, the problem already appears at lower orders. In our calculation to 
order 30, we had to tabulate these matrix elements and compute them independently with 
a C code and the GNU Multiple Precision Arithmetic Library (GMP) which allows arbitrary 
precision floating point operations. The values of the matrix elements were computed 
up to ripjiit = 150 and performing intermediate calculations with 300 significant decimal 
digits. In summary, it turns out that this method is less efficient than the one explained in 
section 0. Nevertheless, apart from serving as a check of the results, we found interesting 
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to explore this quantum mechanical proeedure, since it seems directly generalisable to the 
construction of self-dual configurations in the four-dimensional torus (13), and all the main 
intermediate objects (as (n, s'\U{ki,k2)\p, s)) appear there as well. This is currently under 
study. 
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